Consignes

Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :

La bioinfo c'est : <code>MERVEILLEUX</code>.

N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.

Introduction

Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :

Analyses

Organisation de votre espace de travail

# Je cree differents repertoires afin de stocker les données au fur et à mesure des etapes
mkdir -p ~/Module4-5/EvaluationM4M5-main/DATA/FASTQ  
mkdir -p ~/Module4-5/EvaluationM4M5-main/DATA/CLEANING
mkdir -p ~/Module4-5/EvaluationM4M5-main/DATA/MAPPING
mkdir -p ~/Module4-5/EvaluationM4M5-main/DATA/QC
# Je me deplace dans mon repertoire et valide la presence des dossiers
cd ~/Module4-5/EvaluationM4M5-main
tree ~/Module4-5/EvaluationM4M5-main

Téléchargement des données brutes

Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]

#Tout d'abord reserver les ressources de calcul avec salloc
salloc --cpus-per-task=6 --mem=1G
#charger sra-tools et identifier la version
module avail sra
##sra-tools/2.10.0  sra-tools/2.10.3 
module load sra-tools/2.10.3
#avec fasterq-dump recuperer les données puis les compresser
srun --cpus-per-task=6 fasterq-dump --split-files -p SRR10390685 --outdir ./DATA/FASTQ
srun gzip *.fastq

Verifier les fichiers telechargés en lisant les 6 premieres lignes.

cd ~/DATA/FASTQ
zcat SRR10390685_1.fastq.gz | head -n 6
zcat SRR10390685_2.fastq.gz | head -n 6

Combien de reads sont présents dans les fichiers R1 et R2 ?

module load bc/1.07.1
echo $(zcat SRR10390685_1.fastq.gz | wc -l) /4| bc
echo $(zcat SRR10390685_2.fastq.gz | wc -l) /4| bc

Les fichiers FASTQ contiennent 7066055 reads.

Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

wget ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz

Quelle est la taille de ce génome ?

#Utiliser l'outil seqkit
module load seqkit/0.14.0
seqkit stats GCF_000009045.1_ASM904v1_genomic.fna.gz

La taille de ce génome est de 4215606 paires de bases.

Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

wget ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz

Combien de gènes sont connus pour ce génome ?

zgrep ID=gene GCF_000009045.1_ASM904v1_genomic.gff.gz | wc

4536 gènes sont recensés dans le fichier d’annotation.

Contrôle qualité

J’aimerai avoir quelques informations sur ces 2 fichiers fastq

module load seqkit
srun seqkit stats --threads 1 *.fastq.gz
file format type num_seqs sum_len min_len avg_len max_len
SRR10390685_1.fastq.gz FASTQ DNA 7,066,055 1,056,334,498 35 149.5 151
SRR10390685_2.fastq.gz FASTQ DNA 7,066,055 1,062,807,718 130 150.4 151

On voit que dans les 2 fichiers on a le meme nombre de sequences. En revanche la taille minimale dans le #1 est de 35 alors que dans le #2 elle est de 130

Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit

cd ~/Module4-5/EvaluationM4M5-main/DATA
salloc --cpus-per-task=8 --mem=1G
module load fastqc/0.11.9
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_1.fastq.gz -o QC/ -t 8
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_2.fastq.gz -o QC/ -t 8

La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?

  • Oui
  • Non

car le Phred score median est toujours superieur à 30 comme le montre le graphique Per base sequence quality dans les FastQC Report

Lien vers le rapports MULTIQC

Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?

  • Oui
  • Non

car les reads entre 1_fastq et le 2_fastq n’ont pas la meme longueur minimale (35 versus 130). Il n’y a plus d’adaptateur sur les reads de 1_fastq, mais il reste des Illumina Universal Adapter sur les reads de 2_fastq

Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?

#La taille du génome de reference est de 4215606 paires de bases
#Le nombre de reads est de 7066055 *2
#P= (7066055*2)*150/4215606
#502X

La profondeur de séquençage est de : 500 X.

Nettoyage des reads

Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.

module load fastp
#noter la version de fastp: fastp --version #fastp 0.20.0
srun --cpus-per-task 8 fastp --in1 FASTQ/SRR10390685_1.fastq.gz --in2 FASTQ/SRR10390685_2.fastq.gz --out1 CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz --out2 CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz --html CLEANING/fastp.html --thread 8 --n_base_limit 5 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail  --overrepresentation_analysis --json CLEANING/fastp.json

Les paramètres suivants ont été choisis :

Parametre Valeur Explication
n_base_limit 5 pour eliminer les sequences polyN surrepresentees dans le 1_fastq
cut_mean_quality 30 la mediane des Phred score est superieure à 30. Les reads sont de bonne qualité
cut_window_size 8 fenetre glissante de 8 bases pour le calcul moyen de la qualité
length_required 100 la longueur minimale des reads doit etre de 100
cut tail parametre pour couper ici la queue
overrepresentation_analysis pour eliminer les sequences surrepresentees

Ces paramètres ont permis de conserver 6777048 reads pairés, soit une perte de 4.09% des reads bruts.

Alignement des reads sur le génome de référence

Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].

cd ./MAPPING/
module load bwa/0.7.17
#il faut tout d'abord indexer le genome de reference
srun bwa index GCF_000009045.1_ASM904v1_genomic.fna.gz #GCF_000009045.1_ASM904v1_genomic.fna.fai

#reserver un job allocation adapté
salloc --cpus-per-task=32 --mem=8G

#debuter le mapping à partir des reads filtres
module load bwa mem2/2.2.1
srun --cpus-per-task=32 bwa mem GCF_000009045.1_ASM904v1_genomic.fna.gz ../CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz ../CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz -t 32 > SRR10390685_on_GCF_000009045.1_ASM904v1.sam
#creer le fichier bam
module load samtools/1.10
srun --cpus-per-task=8 samtools view --threads 8 SRR10390685_on_GCF_000009045.1_ASM904v1.sam  -b > SRR10390685_on_GCF_000009045.1_ASM904v1.bam

#trier le fichier bam
srun samtools sort SRR10390685_on_GCF_000009045.1_ASM904v1.bam -o SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam
#indexer le fichier bam trié
srun samtools index SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam


#Enfin je genere un multiqc global
cd ~/Module4-5/EvaluationM4M5-main/DATA/
module load multiqc
srun multiqc -d . -o .

Combien de reads ne sont pas mappés ?

cd ./MAPPING
srun samtools flagstat 
SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam.flagstat

srun samtools idxstats SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam > SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam.idxstats

13571369 reads totaux - 12826829 reads mappés=744540

744540 reads ne sont pas mappés.

Croisement de données

Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:

cd ..
#Chercher le gene trmNF dans le fichier d'annotation du genome complet
zgrep trmNF GCF_000009045.1_ASM904v1_genomic.gff.gz > trmNF.gff3



#Croiser le fichier d'alignement de tous les reads avec la sequence de trmNF et ne recuperer que les reads avec 50% d'overlap
module load bedtools
srun bedtools intersect -a ./MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1.sort.bam  -b trmNF.gff3  -f 0.50 > trmNF_reads.bam

# Tri et indexage des reads
srun samtools sort trmNF_reads.bam -o trmNF_reads.sort.bam
srun samtools index trmNF_reads.sort.bam

srun samtools flagstat trmNF_reads.bam > trmNF_reads.bam.flagstat

2801 reads chevauchent le gène d’intérêt.

Visualisation

Utilisez IGV [7] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.

J’ai eu du mal à generer un fichier indexé .fai de mon genome de reference. Apres plusieurs recherches, orientées par les messages d’erreur je l’ai généré en faisant ceci:

cd ./MAPPING/
zcat GCF_000009045.1_ASM904v1_genomic.fna.gz | bgzip -c > Genome.fa.gz
samtools faidx Genome.fa.gz

Ma capture d’ecran d’IGV:

Attention vue partielle des reads mappés

ImageFinale

References

1. toolkit NS. NCBI SRA toolkit. NCBI, GitHub repository. 2019.
2. Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
3. Zhou Y, Chen Y, Chen S, Gu J. Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90. doi:10.1093/bioinformatics/bty560.
4. Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:13033997. 2013.
5. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.
6. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
7. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative genomics viewer (IGV): High-performance genomics data visualization and exploration. Briefings in bioinformatics. 2013;14:178–92.
 

A work by Migale Bioinformatics Facility

https://migale.inrae.fr

Our two affiliations to cite us:

Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France

Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France